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We show that the problem of finding the primary and secondary characteristic directions of a 
linear lossless optical element can be formulated in terms of an eigenvalue problem related to the 
unimodular factor of the transfer matrix of the optical device. This formulation makes any ac- 
tual computation of the characteristic directions amenable to pre-implemented numerical routines, 
thereby facilitating the decomposition of the transfer matrix into equivalent linear retarders and 
rotators according to the related Poincare equivalence theorem. We explain in detail how this issue 
, arises in the context of stress analysis based on integrated photoelasticity or hybrid methods com- 

D ■ bining photoelastic measurements with analytical stress models and/or numerical Finite-Element 

pL^ ' computations for the stress tensor field. Furthermore we show how our results can be applied when 

algorithms for the reconstruction of the dielectric tensor in the interior of a photoelastic model 
(dielectric tensor imaging) are tested for their stability against noise in the measurement data. For 
the sake of completeness we provide a brief derivation of the basic equations governing integrated 
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photoelasticity. 
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Passive linear optical elements are ubiquitous in the study of interactions between matter and polarized classical [1— 
5] or quantum light [6]. While polarizers attenuate one of two distinguished orthogonal polarization forms, linear 
retarders and rotators alter the state of polarization while preserving the flow of light energy through the device. The 
t> ■ latter two belong to the class of non-absorbing linear elements, where the term "linear" refers to the fact that the 
action of such a device on the state of polarization can be conveniently described by a unitary two-by-two transfer 
matrix. In contrast, a polarizer would be represented by a Hermitean matrix. This description of linear optical 
elements in terms of Hermitean and unitary matrices is called the Jones formalism [7-15]. 

In his treatise on classical light [16], Poincare found that any non-absorbing passive linear optical element could be 
decomposed into basic linear retarders and rotators. By linear retarder we mean a homogeneously anisotropic device, 
such as a piece of appropriately cut crystal, which possesses two preferred axes, called the "fast" and "slow" axis, 
which are perpendicular to each other and differ in the phase velocity of component waves linearly polarized along 
the distinguished directions; as a consequence, light possessing a general elliptic polarization state will accumulate 
a relative phase retardation between the two components and thus change its polarization form. A rotator, on the 
other hand, changes the plane of linearly polarized light by a specified angle; it can be shown easily that this effect 
is due to a phase retardation between the two orthogonal components of circular polarization. Accordingly, a linear 
Oh- retarder is determined by specification of, e.g., the angle of the fast axis, and the relative phase retardation; while a 
single rotation angle is sufficient to specify a rotator. This decomposition of a general non-absorbing optical element 
into retarders and rotators is called the Poincare equivalence theorem (see Ref. [17] for a recent account). 

Such a decomposition proves very useful whenever the determination of the internal stress tensor of a transparent 
medium by means of non-destructive methods is desired. This issue arises in the following context: material scientists 
and engineers wish to gain insight into the stress and strain fields in the interior of a loaded specimen. Amongst the 
many different approaches, three methods are most often used: (1) an educated guess at the analytic mathematical 
form of the internal stress or strain tensor is made; (2) the external load acting on the specimen is systematically 
approximated by forces applying on finitely many points on the surface of the object; and the resulting stress in the 
interior is then computed using numerical schemes generically called Finite Element Methods (FEM) [18-20]; (3) the 
phenomenon of photoelasticity is utilized to gain insight into the internal stress distribution. The term photoelasticity 
refers to the fact that some transparent materials, typically resins or glasses, become birefringent when under external 
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load, while being optically isotropic in the unloaded state. Within certain limits, the relation between the dielectric 
tensor and the stress tensor is linear, and is generically called a stress-optical law. Its typical form has been 
given long ago by Maxwell [21], 

€ij = e Sij + Ci (Jij + C 2 tr a 6ij , (1) 

where Ci, C2 are called stress-optical constants. This one-to-one relation between dielectric tensor and stress tensor 
suggests that one builds a model of an industrial component from an appropriate material (today, plastics are usually 
used) and loads the model just as the real object. The model is then subjected to heat, so as to loosen, and rearrange, 
the molecular bonds in the material; and subsequently cooled down, or "stress-frozen" , upon which the stress pattern 
remains locked inside the material [22-24]. To determine these patterns, destructive and non-destructive methods are 
available. 

A particular example of non-destructive evaluation has been termed "Integrated Photoelasticity" [25]. Here, polar- 
ized light is sent through the specimen at many different angles, and the change in polarization form is registered for 
each light ray (pixelwise), utilizing appropriate combinations of polariscopes and digital cameras. To each ray passing 
through a specified point, and along a specified direction, a unitary transfer matrix U can be assigned, which describes 
how the polarization form changes along the given ray. In this way we obtain a map from the set of all (reasonably 
smooth) dielectric tensor, or stress tensor, fields in the interior of the specimen, to the collection of transfer matrices, 
gathered for all necessary points and directions. In the theory of Inverse Problems and Mathematical Tomography 
[26-31], this map is commonly called the forward problem. The associated inverse problem consists of reconstructing 
the interior dielectric tensor, or stress tensor, from a given collection of transfer matrices. 

The transfer matrices U are not observable directly. Rather, they are known functions of the global phase of the light 
beam, and three so-called "characteristic parameters" [25], two of which may be regarded as polarization directions, 
while the third one has the meaning of a "characteristic phase retardation". Their operational meaning is as follows: 
the primary characteristic directions determine those planes of linear polarization at the entry into the medium for 
which the state of polarization of the emergent beam is again linear. The secondary characteristic directions determine 
the planes of linear polarization of the emerging light, if the incident light was linearly polarized in the primary 
characteristic directions; in general, they differ from the primary ones. There are always two orthogonal primary and 
two orthogonal secondary characteristic directions such that light which is linearly polarized along the two primary 
directions travels with different phase velocities. As a consequence, both waves emerge with a phase difference — the 
characteristic phase retardation. The characteristic parameters [25] of a linear lossless device are usually taken to be 
simple linear combinations of the angles specifying one of the primary characteristic directions, and the angle of the 
associated secondary characteristic direction; and the characteristic phase retardation (see section IV). 

Traditionally, the forward problem in integrated photoelasticity is formulated directly in terms of the bulk stress 
tensor; however, more recent attempts [32] have shown that it may be advantageous to first determine the dielectric 
tensor inside the model by tomographic means, after which the stress tensor can be computed via the stress-optical 
law in cq. (1). The solution to the inverse problem for the two-dimensional (2D) case, i.e., when the specimen is a 
thin slab, and the stress tensor inside possesses only two principal directions, is known long since [3, 22-25, 33-36]. 
The general solution to the three-dimensional (3D) inverse problem is not known to date, although it has been shown 
recently [32] that, in the limit of weak optical anisotropy, the 3D inverse problem can be mathematically reduced to 
six independent 2D inverse problems. 

Being one of the oldest methods of experimental stress analysis, photoelasticity has been somewhat overshadowed 
by FEM methods over the past two or three decades, but has seen a recent revival with applications in silicon wafer 
stress analysis, rapid prototyping, fiber optic sensor development, and image processing [37]. As photoelasticity can 
be regarded as a natural complement to FEM numerics, hybrid methods attempting at combining both approaches 
are becoming popular. It then becomes a viable task to compare FEM results with those obtained by photoelastic 
tomography. 

These considerations describe the context in which the method of determining the characteristic directions of a 
linear optical element, as given in this paper, is expected to prove useful. Specifically, there are four applications of 
our results in integrated photoelasticity and associated hybrid methods: 

1. Testing the validity of analytic stress models: Here, we start with an educated guess about the analytic form of 
the bulk stress tensor field; then compute the associated dielectric tensor via the stress-optical law in eq. (1); and 
use this to compute the transfer matrices of the forward problem numerically, using eq. (12) in section II below. In 
order to facilitate a comparison with the characteristic parameters obtained from a direct photoelastic evaluation 
of a stress-frozen model, a decomposition of the numerically obtained transfer matrices into equivalent retarders 
and rotators must be performed for each light ray. Such a decomposition has to evaluate the characteristic 
directions of a given transfer matrix U first; subsequently, the associated characteristic phase retardation follows 
automatically, as will be shown in section III. It is here where our method of determining the characteristic 
directions is needed. 
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2. Testing the validity of a FEM-based calculation: instead of an analytic stress model we might want to compare a 
bulk stress tensor field evaluated numerically using FEM, with photoelastic data. As in the last item, the stress 
tensor leads to a dielectric tensor and subsequently to a collection of numerically computed transfer matrices 
via the forward problem, cq. (12). Again, the comparison with actual measurements made on a photoelastic 
specimen requires the decomposition of these transfer matrices according to the equivalent optical model. 

3. Iterative solutions of the inverse problem in integrated photoelasticity: the integral equation (12) determining 
the transfer matrices in terms of the dielectric tensor field is non-linear; this accounts for the fact that the full 
solution to the general 3D inverse problem of determining the dielectric tensor in terms of the transfer matrices 
is not yet known. Iterative numerical schemes to solve cq. (12) for may be conceived. At each stage of 
iteration, the intermediate result for the dielectric tensor can be used to compute a collection of associated 
transfer matrices, whose characteristic parameters, in turn, may be compared with actual photoelastic data in 
order to check the convergence of the iterative algorithm. At this last stage we again need the decomposition 
of a given transfer matrix in terms of the equivalent optical model. 

4. Stability of inversion algorithms on noisy data: the reliability of a reconstruction depends on how stable the 
algorithm is with respect to finite errors in measurement data. This stability can be tested on artificial analytic 
stress models, as follows: the stress model can be converted into a collection of transfer matrices as explained 
in item 1.) above. We then decompose each U into characteristic parameters ; but instead of feeding these 
data directly into the inversion we subsequently add some numerical noise to the characteristic parameters, thus 
simulating errors in the measuring apparatus. These modified data are then fed into the proposed algorithm to 
check how far the new reconstruction deviates from the one without noise in the data. A visual example of this 
procedure is given in Fig. 1 for a standard reconstruction algorithm (filtered back-projection [30]). 




(a) (b) 



FIG. 1: A cylinder is subject to axial load; the resulting deformation is modelled as a simulated displacement field on the 
surface and in the interior. Hooke's law gives the relation between the resulting strain and stress tensors; from the stress- 
optical law, eq. (1), the dielectric tensor is then obtained. This is an example of an artificial stress model giving rise to a 
simulated dielectric tensor which can be used to test photoelastic reconstruction algorithms as follows: a plane intersecting the 
cylinder at an angle of 22° relative to the symmetry axis of the cylinder is chosen in the figure. The unit normal vector to this 
plane is rj. The component of the dielectric tensor in the direction of rj is e vv — e(r],ri). In subfigure (a), its relative deviation 
(e?7?7 — e)/e = A m from the scalar isotropic value e is plotted. In subfigure (b), the simulated dielectric tensor £»j is used to 
compute a collection of transfer matrices U via the forward problem, eq. (12). The transfer matrices are then decomposed 
according to the equivalent optical model, using the method introduced in section III; subsequently, Gaussian noise is imposed 
on the characteristic parameters, eq. (37), thus simulating measurement errors. The noisy parameters then recombine, again 
via the equivalent optical model, to give new, "noisy", transfer matrices, from which the "noisy" dielectric tensor can be 
reconstructed via a filtered back-projection algorithm [30, 31]. The result of this reconstruction, for 60 x 60 pixels and 90 
scans around the object, is plotted in subfigure (b). It is seen that the noise in the characteristic parameters, and the noise in 
the global phase (j>, eq. (16), give rise to a halo- like pattern around the object, in addition to the speckled appearance of the 
reconstruction. - The comparison of Figures (a) and (b) provides a visual example of how our decomposition method can be 
used in order to test the stability of a reconstruction algorithm (here: filtered back-projection) against measurement errors. 

In Ref. [17] we have given a detailed account of the Poincare equivalence theorem and its application to construct 
the equivalent optical model. In that paper, the decomposition of a given transfer matrix U into retarders and rotators 
was accomplished in a somewhat pedestrian fashion, by parameterising the (three-dimensional) manifold of SU (2) 
matrices in a suitable way, and then applying elementary trigonometric methods. While mathematically correct, the 
parameterisation of the SU(2) manifold as used in Ref. [17] is inconvenient for computer-based algorithms in that 
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subsets of the boundary of the parameter space of this manifold must be identified in order to obtain a one-to-one 
relation between parameters and matrices. This non-trivial topology introduces complications when used in the 
numerical schemes outlined above. In the present paper we accomplish the same decomposition in a more elegant 
way: we show how to reformulate the problem of finding the "characteristic" data specifying the equivalent optical 
model for a given transfer matrix in terms of an eigenvalue problem associated with the unimodular factor of this 
matrix. The new method is therefore more suitable for an actual numerical determination of the characteristic data, 
as we can immediately make use of pre-implemented numerical routines for eigenvalue problems. 

The plan of the paper is as follows: in section II we briefly discuss the derivation, and physical content, of the basic 
equations governing integrated photoelasticity, leading to the integral equation (12) for the transfer matrices in terms 
of the dielectric tensor, which can be regarded as the basic equation governing integrated photoelasticity. In section III 
we introduce our new scheme for determining characteristic parameters for a given optical transfer matrix in terms of 
an eigenvalue problem. In section IV we explain the relation of the results so obtained with the Poincare Equivalence 
Theorem, which was discussed more thoroughly in Ref. [17]; and we define the equivalent optical model, and the 
associated characteristic parameters, in terms of the characteristic directions determined by the method introduced 
in section III. In section V we summarize the paper. 



II. THE EQUATIONS GOVERNING INTEGRATED PHOTOELASTICITY 

We consider a non-magnetic material which is transparent for optical wavelengths; and which is optically isotropic 
(ey = eSij) when unloaded but becomes birefringent when under external load. Plastics and certain resins exhibit 
these properties. 

As explained in Ref. [32], the spatial variation of the optical inhomogeneity in the material — acquired under 
loading — can be characterised by the trace (1/3) tre of the dielectric tensor; a heuristic length scale lo then gives 
the distance over which the relative change in (1/3) tre is comparable to one. Furthermore, since the loaded medium 
is optically anisotropic at each point, two preferred polarization directions [1-3, 38, 39] exist for each direction of prop- 
agation. If the medium were homogeneously anisotropic, like a crystal, these preferred directions would be constant 
through the whole medium; in photoelasticity, however, the anisotropy in will vary at each point in the material, so 
that another heuristic scale l p , denoting the distance along which the relative change of a given preferred polarization 
direction is comparable to one, may be defined. If these two scales are substantially larger than the wavelength of 
the light passing through the material, Iq,I p 3> A, the propagation of light through the specimen can be described in 
the geometrical- optical approximation [1, 39], where the complex electric field of the light beam may be written in the 
form 

E(x,t) = E(x)e i * (x)_<wt . (2) 

Here, the eikonal <f)(x) describes a locally-plane wave, with local wave vector V0(x), and local amplitude E(x). By 
assumption, both of these quantities vary weakly on the length scale of a wavelength A. It is then convenient to 
introduce a dimensionless scale [40], 

a = maX {^} ' (3) 

in terms of which the geometrical-optical limit is simply characterised by a < 1. 

Furthermore, an explicit measure of anisotropy may be given by the relative deviation of e from its isotropic 
value [32], 

Aij = e JizJiil , (4a) 

/3 ~ max || Aij || , (4b) 

where the global maximum j3 characterises the magnitude of anisotropy in the dielectric tensor. In principle, this 
inhomogeneous anisotropy would lead to a continuous splitting of light rays in the medium [40], due to the fact 
that two distinct phase velocities, and two distinct ray velocities, for any given propagation direction exist at each 
point [1, 2, 39]. However, experimentally, no ray splitting is seen in typical photoclastic measurements; rather, light 
rays propagate, for all practical purposes, along straight lines through the object, and the only effect of the optical 
tensors on the light beam is to rotate the preferred polarization directions. This observational fact indicates that, in 
the context of (industrially relevant) photoelasticity, the anisotropy of is so weak as to permit a replacement of 
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the two actual rays by one single, "effective" ray which is obtained from the isotropic part (1/3) tre of the dielectric 
tensor alone. The mathematical condition for this to be true has been given in Rcf. [40], 







a 



< 1 



(5) 



In this case, the local wave vector V</> in eq. (2) may be taken as globally constant, V<p = k, so that (2) describes a 
strictly plane wave 

E(x) e *-*-*"* . 



E(x, t) 



(6) 



Only the weakly varying amplitude E(x) then encodes the structural inhomogencities in the material. 

In the limit of the geometrical-optical approximation a«l, and under the condition of negligible ray splitting (5), 
the condition of weak anisotropy {3 <C 1 is automatically satisfied. In this case the electric field may be taken as 
transverse to the propagation direction k = k/k of the light beam [32]. Upon inserting the trial solution (6) into the 
wave equation 



AE - fi eE = 



(7) 



and retaining only first-order spatial derivatives — corresponding to the geometrical-optical limit — we arrive at an 
equation of the form 



k x x e) - ^-|v x(kxE) + kx(VxE)|+ Mo" 2 eE = . 



(8) 



Here, u = uj/k is the phase velocity in the unloaded material, and e is the dielectric tensor. The longitudinal 
component of eq. (8), obtained by projection onto the unit vector k, can be neglected in the geometrical-optical limit. 
In a coordinate system for which k is along the z axis, eq. (8) becomes 

(9) 



(10) 
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= l x 


A21 A22 




E2 



with Aij as defined in eq. (4a). The solution of (9) can be expressed in terms of a transfer matrix U such that 



Ei(z) 
E 2 {z) 



= U(z,z ) 



Ei(z ) 
E 2 (z ) 



where U satisfies 



d 7T 

— U(z,z ) =i-Aj_(z)U(z,z ) 



1 
1 



(11) 



U(z ,z ) 

T n < system (11) is equivalent to an integral equation 

7T 

U(z,z ) = ± 2 + ij I dzi A±(z 1 )U(z 1 ,z ) 



(12) 



where A± denotes the matrix of transverse components of A^ as they appear in eq. (9), and A is the wavelength in 
the unloaded material. Eq. (12) is the basic equation governing integrated photoelasticity: it determines the transfer 
matrices U for a given light ray, as functions of the anisotropic part A^ of the dielectric tensor in the interior of the 
medium. It can be shown that U must be unitary, preserving the norm of the complex electric field vector. This is 
clear on physical grounds, as preservation of the norm of E implies preservation of intensity, so unitarity here just 
expresses energy conservation of the light beam. This is just to be expected, since we have assumed a non-absorbing 
medium. 

Eq. (12) is manifestly non-linear in U, since the transfer matrices appear under the integral on the right-hand 
side. This can be seen by iteratively inserting the left-hand side of eq. (12) into the right-hand side, leading to the 
Born-Neumann series 

Zl 



U(z,z ) = 1 2 + (ij^j J dz x A ± ( Zl ) + (i^j J dz x A ± ( Zl ) J dz 2 A ± (z 2 ) + 



(13) 



This nonlinearity provides one of the major mathematical challenges in the theory of integrated photoelasticity. As a 
consequence, the inverse problem of determining A^ in terms of a collection of transfer matrices is highly non-trivial. 
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III. CHARACTERISTIC DIRECTIONS OF LINEAR OPTICAL ELEMENTS 

In the previous section we have seen that, in the theory of integrated photoelasticity, a transfer matrix U can be 
assigned to each light ray which passes a photoelastic material through a given point, and along a specified direction. 
This means that the associated light path may be regarded as a linear lossless optical element which acts on the given 
light ray by changing its polarization form. As explained in the introduction, there are several occasions in photoelastic 
stress analysis, hybrid FEM-photoelastic methods, and methods determining the numerical stability of numerical 
algorithms aiming at reconstructing the dielectric tensor from the collection of — tomographically determined — transfer 
matrices, where the decomposition of a given transfer matrix U according to the equivalent optical model is an issue. 
As will be shown in this section, the major step in this task is to determine the characteristic directions of the transfer 
matrix U. 

As shown in [17], a lossless linear passive optical device possesses in general two so-called primary characteristic 
directions [25] w m — (cos-f m ,smj m ),m = 1,2, in the plane perpendicular to the entry of the optical element, which 
have the following significance: if a light beam at the entry is plane-polarized in one of the directions w m (our 
convention is such as to define the direction of polarization along the electric displacement field D) it will leave the 
device again in a state of plane polarization, with the plane oriented along unit vectors w' m — (cos j' m , sin j' m ), m = 1,2, 
called the secondary characteristic directions. In contrast, for any direction other than the beam at exit will in 
general be elliptically polarized. The two primary as well as the two secondary characteristic directions are always 
perpendicular to each other, so that it suffices to specify the angle 7 and 7' of the first elements w\ — (cos 7, sin 7) and 
w[ = (cos 7', sin 7'), respectively; the second element W2,w' 2 is then determined up to a sign. Since the polarization 
state at the exit is again linear the optical device, represented by the unitary matrix U, must act on the real polarization 
vector w m according to 

Uw m = e i *™w' m , (14) 

where $1, $2 are the phases picked up by the light beam entering along w\,W2, respectively. Our goal is to determine 
a consistent choice of primary and secondary characteristic vectors, together with appropriate values for the phases 
$ m , from a given transfer matrix U. 

Since U is unitary, its determinant is a unimodular number 

exp(2i0) = det U , (15) 

hence we can factorize U into 

U = exp(#) S , (16) 

where S is now a unimodular unitary matrix, det S = 1. The choice of S is not unique, since both S and — S satisfy 
det S — 1. The phase <f> can be computed from (15) modulo tt, the ambiguity in sign obviously related to the double- 
valuedness ±5 of the SU (2) factor. We therefore need to stipulate an explicit convention for the two possibilities in 
the factorization: we choose <j) to be the smallest possible non-negative solution of (15). Then S is uniquely determined 
by eq. (16). 

We can now rewrite (14) as 

Sw m = e i *'™w' m , $:„ = $ m -0 . (17) 

In principle we can determine the angles 7 and 7' of the primary and secondary directions w m and w' m by parametrising 
the manifold of SU (2) matrices S in a suitable way and then using elementary trigonometric relations to express these 
angles in terms of the coordinates on the SU{2) manifold, as was done in [17]. However, a method that docs not 
require an explicit coordinate chart on the SU (2) manifold will be presented now: 

Suppose that eq. (14) holds. Then (17) holds as well, and on taking the complex conjugate of the latter equation 
we obtain 

S*w m = e- i ^w' m . (18) 

On eliminating w' m from eqs. (17) and (18) we find that the directions w m are real eigenvectors of S T S, 

(S T S)w m = e 2 ^'w m , (19) 

where the superscript T denotes a matrix transpose. We therefore need a method to obtain real eigenvectors from a 
complex matrix of the form S T S, where S is an element of SU(2). To this end we show that S T S commutes with its 
complex conjugate (S T S)*: any SU(2) matrix S can be represented in the well-known form 



(20) 
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so that a short computation gives 



(S T S) {S T S)* = {4 3m 2 (ab) + \a 2 + b* 2 \ 2 } 1 2 



(21) 



Since the right-hand side is real it follows that the left-hand side is equal to its complex-conjugate; as a consequence, 
the commutator 



(s T s), (s T sy 







(22) 



must vanish. 

This relation is significant, since it can be used, in turn, as a starting point to determine the characteristic directions 
in an elegant way: given any linear lossless device with unitary Jones matrix U = cxp(i</>) S, the SU(2) factor 5* will 
satisfy relation (22). This relation implies that the commutator 



9te5 T 5, 3mS T S 



l{s T s + (s T sy}, l l {s T s-(s T sy} 



(23) 



-±,[s T s,(s T sy] =0 



must vanish as well. The real and imaginary parts of S T S are symmetric, since S T S is. As a consequence of the 
commutativity (23), both of these matrices share the same (orthogonal) system of eigenvectors w m which must be 
real since the real and imaginary parts are so, 



me (S T S) w v 



3m (S T S) w m = j m w m 



(24) 



It follows that w m are real eigenvectors of S T S as well, with eigenvalues r m + ij m . On the other hand, since S T S is 
unitary, its eigenvalues must be the unimodular numbers exp(2i$' m ) appearing in eq. (19), so that 



exp(2i$' m ) = r m + ij„ 



(25) 



The result (24) shows that the characteristic directions w m can be obtained as the real eigenvectors of the matrices 
£Re S T S or 3m S T S describing the optical element. Its significance lies in the fact that the process of finding the 
characteristic directions from a given transfer matrix U becomes amenable to well-established numerical routines 
for general eigenvalue problems. This problem arises e.g. in integrated photoelasticity, or more generally, in any 
effort to reconstruct the dielectric tensor inside a transparent but inhomogencously anisotropic optical device from 
tomographic measurements [32]. - This outlines the principle of our method to obtain the characteristic directions 
of any transfer matrix U. To finish our discussion we now show how to fix the ambiguity in signs of the eigenvectors 
appearing in (19) and (24), and the associated phase ambiguity, in a consistent way: 

The numerical routine will deliver two eigenvectors w m , but there are four possible choices 



(wi, w 2 ) 



(wi,-w 2 ) , {-wi,w 2 ) , (-101,-102) 



(26) 



for the signs. We thus need to agree on a convention to pick a system from (26): we first choose from ±w\ the vector 
which makes an angle with the x axis whose modulus does not exceed 7r/2, so that the ^-component of this vector is 
always non-negative. Without loss of generality we may assume this to be true for +w\. Next we choose from ±W2 
that vector which makes the system (101, ±102, 63) right-handed, where it is assumed that the light beam propagates 
along e 3 towards positive z- values. Without loss of generality we may assume that this is satisfied by +w 2 - The 
phases $>' m can now be obtained from (25), but obviously only modulo ir, 



(27) 



Accordingly we can determine the associated secondary characteristic directions w' m from (17), but only up to a sign, 
due to the phase ambiguity. We therefore have four possibilities 



(w[,w' 2 ) , (w[,~w 2 ) , (-w[,w 2 ) , {-w'^-w'2) 



(28) 



for the secondary system. We now impose two conditions similar to those that made the choice of w m unique: firstly, 
we require that the suitable candidate ±w[ for the first element makes an angle with w\ whose modulus is not larger 
than 7r/2. Assuming that this is the case for w[, we must have 



w[ ■ w\ > , I7' — 7I < 



(29) 
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where w[ — (cos 7', sin 7'). We then still have the ambiguity of ±w' 2 ; our second condition now is to select this sign 
so that the vector triad (w^, ±w' 2 , e^) is right-handed; we may assume that +w' 2 is the correct choice. 

We have now fixed the ambiguous signs of the characteristic directions; as a consequence, the phases <fr' m are 
determined by eq. (17) up to multiples of 2tt. This last indeterminacy is intrinsic and cannot possibly removed. Thus 
we stipulate to let the phases §' m take values in the interval — n < <fr' m < tt, it being understood that the values of 
and $2 are different; for, if they were equal, they would have been part of the phase <j) which was extracted out of U 
in eq. (15). As a consequence, 

-2tt < <5>[ + $ 2 < 2tt . (30) 

S can now be represented in terms of primary and secondary characteristic directions, and associated phases, as 

S = K) exp(i$i) ( Wl \ + K) exp(i$ 2 ) (w 2 \ , (31) 

where we have denoted (column) vectors as \w) and (row) covectors as (w\, reminiscent to quantum-mechanical 
conventions. 



IV. RELATION TO THE POINCARE EQUIVALENCE THEOREM 

Finally we show how the representation (31) of S in terms of characteristic directions and phases is related to the 
decomposition of a lossless linear optical element according to the Poincare equivalence theorem [16]: to this end we 
represent eq. (31) in the basis ei,e2 of Cartesian coordinate vectors pertaining to the laboratory frame: using the 
notation {ei\S\ej) = Sij we see that (31) takes the form 

2 

Sij= Yl J Um 2 R(l)m 2j , (32) 

mi, 7712 = 1 

where 

R{~7 )im — ( e -i\ w m) J ^m\m 2 = e m $m 1 m 2 > ^(oOmj = ( w m\ e j) 

i2( 7 ) = 



mim2 

(33) 



cos 7 sin 7 
-sin 7 cos 7 

using the notation conventions of [17]. The vectors in the pairs (wi,w 2 ) and (w^w'2) are orthogonal, and have been 
constructed to make a right-handed system together with e%. It follows that the matrices R("f), R(— 7') are proper 
rotation matries having unit determinant, i.e. elements of SO(2). Then, since det5 = 1 it follows that det J' = 1, 
which implies that the eigenvalues §' m must sum up to a multiple of 27T, <&[ + $ 2 = 2ttN. But, according to (30) this 
restriction can be made stronger, 

+ $ 2 = • (34) 
Finally, on multiplying (32) with exp(i^) we find on using (16) and (34) that 

U = R(-j') J(0, S) i?( 7 ) , (35a) 
J(0,(5) = diag(exp(-i^/2),cxp( J( 5/2)) , -7; = K+4> ■ (35b) 

We recognize that J(0, S) is the Jones matrix of a linear retarder whose fast axis, for 5 > 0, coincides with the x axis of 
the laboratory system, so that light plane-polarized along e\ [e-i) accumulates a relative phase —8/2 {8/2) on passing 
through the device, without changing its linear polarization form, or the orientation of the plane of polarization. On 
using the fact that the transfer matrix of a linear retarder with fast axis making a nonvanishing angle 7 with the x 
axis is given by 

J{ 1 ,5) = R{- 1 )J{Q,5)R{ 1 ) , (36) 
we can rewrite (35a) in the equivalent forms 



U = J( 7 ', 8) R{-i + 7) = R{-i + 7) J( 7 , 8) 



(37) 
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The decompositions (35), (37) express the fact that any linear lossless optical device can be replaced by a sequence 
of one linear retarder and one or two appropriate rotators, at least as far as its optical properties are concerned. The 
fictitious optical device comprised of these retarders and rotators is called the equivalent optical model. The three 
quantities 7,7', 6 are commonly called the characteristic parameters of the equivalent optical model [25], where the 
angle —7' + 7 specifying the equivalent rotator is the same for both forms in eq. (37). In the first (second) form, 
7' (7) is the angle between the customarily selected primary characteristic direction — the fast axis of the equivalent 
retarder — and the x-axis; while 6 is the characteristic phase retardation of the equivalent retarder in both forms. - 
The physical and mathematical content of (35), (37) is called the Poincare equivalence theorem. The decompositions 
as given above coincide with the forms given in [17]. 



V. SUMMARY 



We have presented a method to determine the primary and secondary characteristic directions of a linear lossless 
optical device from an eigenvalue problem formulated in terms of the unimodular factor of the transfer matrix of 
the optical clement. This approach is conceptually more elegant than methods using explicit parametrisations of the 
manifold of SU(2) matrices, and is furthermore amenable to pre-implcmented numerical routines, thus making the 
decomposition of the transfer matrix in terms of equivalent linear retarders and rotators numerically more convenient. 
This important issue arises in the context of stress analysis based on integrated photoelasticity or hybrid methods 
combining photoelastic measurements with analytical stress models and/or numerical Finite-Element (FEM) evalu- 
ations of the stress tensor field. In addition, we have given a visual example of how our results can be used to test 
the stability of reconstruction algorithms for the dielectric tensor in the interior of a photoelastic model on noisy 
measurement data. A brief derivation of the basic equations governing integrated photoelasticity has been provided. 
The relation of our results to the associated Poincare equivalence theorem has been explained. 
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